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Abstract: Assumed stress hybrid methods are known to improve the per- 
formance of standard displacement-based finite elements and are widely 
used in computational mechanics. The methods are based on the Hellinger- 
Reissner variational principle for the displacement and stress variables. This 
work analyzes two existing 4-node hybrid stress quadrilateral elements due 
to Pian and Sumihara [Int. J. Numer. Meth. Engng, 1984] and due to Xie 
and Zhou [Int. J. Numer. Meth. Engng, 2004], which behave robustly in 
numerical benchmark tests. For the finite elements, the isoparametric bilin- 
ear interpolation is used for the displacement approximation, while different 
piecewise-independent 5-parameter modes are employed for the stress ap- 
proximation. We show that the two schemes are free from Poisson-locking, 
in the sense that the error bound in the a priori estimate is independent of the 
relevant Lame constant A. We also establish the equivalence of the methods 
to two assumed enhanced strain schemes. Finally, we derive reliable and ef- 
ficient residual-based a posteriori error estimators for the stress in L^-norm 
and the displacement in if^-norm, and verify the theoretical results by some 
numerical experiments. 
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1. Introduction 



Let C be a bounded open set with boundary T = ToUTn, where meas(r z))>0 
The plane linear elasticity model is given by 

—diver = f in 

(T = Ce{u) in (1.1) 

uIfb = 0,crn|r^ = g, 

where cr E Mgym denotes the symmetric stress tensor field, u G the displacement 
field, e{u) = (Vu + V^u) /2 the strain, f G the body loading density, g G 
the surface traction, n the unit outward vector normal to T, and C the elasticity 
modulus tensor with 

Ce{u) = 2fi£{u) + Adivu I, 
I the 2x2 identity tensor, and yU, A the Lame parameters given hy fi = 2(^u] ' 

A = (i+,gi_2.) for plane strain problems and by /i = X = jYT^friT) 

plane stress problems, with < u < 0.5 the Poisson ratio and E the Young's 
modulus. 

It is well-known that the standard 4-node displacement quadrilateral element 
(i.e. isoparametric bilinear element) yields poor results at coarse meshes for prob- 
lems with bending and suffers from "Poisson locking" for plane strain problems, at 
the nearly incompressible limit (A — )■ oo as — )■ 0.5). We refer to [[II for the mathe- 
matical characteristic of locking. To improve the performance of the isoparametric 
bilinear displacement element while preserving its convenience, various methods 
have been suggested in literature. 

The method of incompatible displacement modes is based on enriching the stan- 
dard displacement modes with internal incompatible displacements. A representa- 
tive incompatible displacement is the so-called Wilson element proposed by Wil- 
son, Taylor, Doherty, and Ghaboussi [29[. It achieves a greater degree of accuracy 
than the isoparametric bilinear element when using coarse meshes. This element 
was subsequently modified by Taylor, Wilson and Beresford [[271 . and the modified 
Wilson element behaves uniformly in the nearly incompressibility. In [[T4|. Le- 
saint analyzed convergence on uniform square meshes for Wilson element. He and 
Zlamal then established convergence for the modified Wilson element on arbitrary 
quadrilateral meshes [15J. In [26J, Shi established a convergence condition for the 
quadrilateral Wilson element. In [[34l . Zhang derived uniform convergence for the 
modified Wilson element on arbitrary quadrilateral meshes. 

The assumed-stress hybrid approach is a kind of mixed method based on the 
Hellinger-Reissner variational principle which includes displacements and stresses. 
The pioneering work in this direction is by Plan [|T6l, where the assumed stress field 
assumed to satisfy the homogenous equilibrium equations pointwise. In [[TTl Plan 
and Chen proposed a new type of the hybrid-method by imposing the stress equi- 
librium equations in a variational sense and by adopting the natural co-ordinate for 
stress approximation. In [ 18[ Plan and Sumihara derived the famous assumed stress 
hybrid element (abbreviated as the PS finite element) through a rational choice of 
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stress terms. Despite of the use of isoparametric bilinear displacement approxima- 
tion, the PS finite element yields uniformly accurate results for all the numerical 
benchmark tests. Plan and Tong [20] discussed the similarity and basic difference 
between the incompatible displacement model and the hybrid stress model. In the 
direction of determining the optimal stress parameters, there have been many other 
research efforts [[11 [30l [3D [32l [3l. In [[30l 123, Xie and Zhou derived robust 
4-node hybrid stress quadrilateral elements by optimizing stress modes with a so- 
called energy-compatibility condition, i.e. the assumed stress terms are orthogonal 
to the enhanced strains caused by Wilson bubble displacements. In [36 [ a conver- 
gence analysis was established for the PS element, but the upper bound in the error 
estimate is not uniform with respect to A. So far there is no uniform error analysis 
with respect to the nearly incompressibility for the assumed stress hybrid methods 
on arbitrary quadrilateral meshes. 

Closely related to the assumed stress method is the enhanced assumed strain 
method (EAS) pioneered by Simo and Rifai [l25l. Its variational basis is the Hu- 
Washizu principle which includes displacements, stresses, and enhanced strains. It 
was shown in [25] that the classical method of incompatible displacement modes is 
a special case of the EAS-method. Yeo and Lee [33] proved that the EAS concept 
in some model situation is equivalent to a Hellinger-Reissner formulation. In [24], 
Reddy and Simo established an a priori error estimate for the EAS method on paral- 
lelogram meshes. Braess [[3l re-examined the sufficient conditions for convergence, 
in particular relating the stability condition to a strengthened Cauchy inequality, 
and elucidating the influence of the Lame constant A. In [4], Braess, Carstensen 
and Reddy established uniform convergence and a posteriori estimates for the EAS 
method on parallelogram meshes. 

The main goal of this work is to establish uniform convergence and a posteriori 
error estimates for two 4-node assumed stress hybrid quadrilateral elements: the PS 
finite element by Plan and Sumihara HTSl and the ECQ4 finite element by Xie and 
Zhou [30]. Equivalence is established between the hybrid finite element schemes 
and two EAS proposed schemes. We also carry out an a posteriori error analysis for 
the hybrid methods. 

The paper is organized as follows. In Section 2 we discuss the uniform stability 
of the weak formulations. Section 3 is devoted to finite element formulations of 
the hybrid elements PS and ECQ4 and their numerical performance investigation. 
We establish the uniformly stability conditions and derive uniform a priori error 
estimates in Section 4. Equivalence between the hybrid schemes and two EAS 
schemes is discussed in Section 5. We devote Section 6 to an analysis of a posteriori 
error estimates for the hybrid methods and verification of the theoretical results by 
numerical tests. 

2. Uniform stability of the weak formulations 

First we introduce some notations. Let L^(T; X) be the space of square integrable 
functions defined on T with values in the finite-dimensional vector space X and with 
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norm being denoted by 1 1 ■ | |o,t- We denote by H^{T] X) the usual Sobolev space 
consisting of functions defined on T, taking values in X, and with all derivatives of 
order up to k square-integrable. The norm on H''(T] X) is denoted by 1 1 • 11^^ : = 
Cl2o<j<k I^Ijt)^^^' wi'^h I ■ \k,T the semi-norm derived from the partial derivatives 
of order equal to k. When there is no conflict, we may abbreviate them to 1 1 ■ | |fc and 
\ ■ \k- Let LK^I) be the space of square integrable functions with zero mean values. 
We denote by Pk{T) the set of polynomials of degree less than or equal to k, by Qk 
the set of polynomials of degree less than or equal to k in each variable. 

For convenience, we use the notation a < b to represent that there exists a 
generic positive constant C, independent of the mesh parameter h and the Lame 
constant A, such that a < Ch. Finally, a ^ 6 abbreviates a <b < a. 

We define two spaces as follows: 

V := H^ny = {uE H\nf : u|r, = 0}, 
L2(^];M2x^), zf meas{Tr,)>0, 



{rGL2(^];M2x^J:^^rrrfx = 0}, T 



N 



where L ^lym) denotes the space of square-integrable symmetric tensors with 



the norm 1 1 • | |o defined by | |t| |g := J^r : rdx, and trr := Th + T22 represents the 
trace of the tensor r. Notice that on the space V, the semi-norm | ■ |i is equivalent 
to the norm 1 1 • 1 1 1 . 

The Hellinger-Reissner variational principle for the model (11.11) reads as: Find 
(cr, u) e S X V with 

a(cr,T) - / T : e(u)rfx = foraUreS, (2.1) 
Jn 

(T : £(v)rfx = F(v) for aU v G V, (2.2) 



n 



where 



n 



'n Jtm 



Here and throughout the paper, cr : t = J2'ij=i ^nd := r — ^trrl. 

The following continuity conditions are immediate: 

a('^)'r) ^ Ikllolkllo, cr, T G S, (2.3) 

/ T : £(v)c/x < ||t||o|v|i, reE,vEV, (2.4) 
Jn 

W <(l|f|l-i + llglLirJ|v|i, vGK (2.5) 
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According to the theory of mixed finite element methods [I6ll3, we need the fol- 
lowing two stability conditions for the well-posedness of the weak problem (12.11) - 
(EJ). 

(Al) Kernel-coercivity: For any t E Z := {t E : J^t : e{^v)dx = for all v E 
V} it holds 

I|t||o < a(T,T); 
(A2) Inf-sup condition: For any v G it holds 

|v|i< sup „ ,, ■ 

O^T€S ||T||o 

The proof of (A1)-(A2) utilizes a lemma of Bramble, Lazarov and Pasciak. 

Lemma 2.1. ([5]) For q E L := < ,2/n( ■ / rh it holds 

I -^ol^^J "^J ^ N — V 

Lgdiwdx 

Ikllo < sup . 

vey |v|i 

The following stability result is given in [HI for the model situation F^r = 0. 
Theorem 2.1. The uniform stability conditions (Al ) and (A2) hold. 
Proof. Firstly we prove (Al). Since 

a(T,T) = / ( -^T^ : + —^—-trTtrT] rfx, 
JnK^l^ 4(/i + A) ; 

we only need to prove ||trT||o < IIt^Ho for any t E Z. 
In fact, for t E Z and any v G V^, it holds 

0=1 T:e(v)c/x= / (^trrl + r^) : e(v)c/x 

Jn 2 

= I -trrdivvfix + / : £:(v)(ix. 
Thus, by Lemma 2.1 we obtain 

P^tIIo < sup ^^^^ — — = sup — — — < 2||t^||o. 

vev |v|i vev |v|i 

This implies (Al). For the proof of (A2), let v G and notice £(v) G S. Then 

k(v)|o< sup ,1 ,1 . 

TGS\{0} llTllo 

Hence (A2) follows from the equivalence between the two norms |£(v)|o and |v|i 
onV. □ 

In view of the continuity conditions, (|2.3I) - (|2.5I) . and the stability conditions, 
(A1)-(A2), we immediately get the well-posedness results: 

Theorem 2.2. Assume that f G V^', g G H^^/'^(Tn). Then the weak problem 
(|ZT|)-(I22]) admits a unique solution (cr, u) G S x V such that 

lkllo + |u|i < ||f||-i + llglLir^- 
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3. Finite element formulations for hybrid methods 



3.1 Geometric properties of quadrilaterals 

In what follows we assume that is a convex polygonal domain. Let be a con- 
ventional quadrilateral mesh of fi. We denote by hx the diameter of a quadrilateral 
K E Th, and denote h := maxKeT^ hx- Let Zi{xi,yi), 1 < z < 4 be the four 
vertices of K, and Tj denotes the sub-triangle of K with vertices Zi_i, Zi and Z^+i 
(the index on Zj is modulo 4). Define 

Pk = niin diameter of circle inscribed in Tj. 

l<i<4 

Throughout the paper, we assume that the partition satisfies the following "shape- 
regularity" hypothesis: There exist a constant g > 2 independent of h such that for 

all Ken, 

hx < QpK- (3.1) 

Remark 3.1. As pointed out in l{34\l . this shape regularity condition is equivalent 
to the following one which has been widely used in literature (e.g. hTTi): there exist 
two constants q' > 2 and < 7 < 1 independent of h such that for all K G Th, 

hx < q'p'kj I cos9]^\ < 7 for 1 < i < 4. 

Here p'^ and denote the maximum diameter of all circles contained in K and 
the angles associated with vertices of K. 

Let K = [—1,1] X [—1, 1] be the reference square with vertices Zj, 1 < i < 4. 
Then exists a unique invertible mapping Fk that maps K onto K with Fk{^, r]) G 
Qj{^, Tj) and = Zi,l <i <4: (Figure 1). Here ^, r/ G [-1, 1] are the local 

isoparametric coordinates. 




This isoparametric bilinear mapping (x, y) = Fk{^, rf) is given by 

4 4 

x = ^x,N,{i,7]), y = ^2/,Ar,(e,r/), (3.2) 

i=l i=X 
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where 



4 



4 



We can rewrite (|3.2I) as 

x = ao + ai^ + a2?7 + ais^r/, y = &o + ^1^ + hv + ^12^^, 

where 



(3.3) 



\ ai2 bi2 j 



(I 1 1 1 \ 

-111-1 

-1-11 1 

VI -11-1/ 



( xi yi\ 

X2 1/2 
X3 2/3 



Remark 3.2. Due to the choice of node order (Figure 1), we always have ai > 
0,62 > 0. 

Remark 3.3. Notice that when K is a parallelogram, we have ai2 = &12 = 0, and 
Fk is reduced to an affine mapping. 



Then the Jacobi matrix of the transformation F^ is 

DFK{^,r]) 





dx \ 






drj 




I dy 


dy 




\ a? 


dr] j 





61 + hi2ri 62 + &i2^ 



and the Jacobian of Fx is 

JK{i. V) = det{DFK) = Jo + M + J2V, 

where 

Jo = 0162 — «2^1, Jl = «1&12 — «12&1, -h = «12&2 — a2&12- 

Denote by the inverse of Fk, then we obtain 



95 95 

dx dy 

drj dr] 

dx dy 



DF-' o FK^i^r,) = [DFk)-' 

1 / 62 + hi2i -a2 - ai2^ 



It holds the following element geometric properties: 
Lemma 3.1. ([34]) For any K G T^, under the hypothesis di.iP . we have 



max Jk{^,v) ,0 2 
min Ji^(^,r7) 2p2^ " 2 ^ 



(3.4) 



ip^ < a? + 6? < i/i^, <al + bl< hi, al^ + fe^^ < ^h^. (3.5) 
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In view of the choice of node order (cf. Figure 1), the shape-regular hypothesis 
(13.11) and the relations (|3.5I) . without loss of generality we assume 

\bi\ < ai, \a2\ < &2- (3.6) 

Together with (13.51) . this leads to 

ai ^ 62 ~ hx, max{a2, h} < O^hx). (3.7) 

Notice also that Lemma lSTTI shows 

Jk^Jo^ hi. (3.8) 



3.2 Hybrid methods PS and ECQ4 



This subsection is devoted to the finite element formulations of the 4-node assumed 
stress hybrid quadrilateral elements PS |tl8J and ECQ4 [|30ll . 

Let S/j C S and Vh CV finite dimensional spaces respectively for stress and 
displacement approximations, then the corresponding finite element scheme for the 
problem (|2.1I)(|2.2I) reads as: Find {crh, u/i) G x Vh, such that 



(3.9) 
(3.10) 



a{(Th, t) — / T : £:(u/i)(ix = for all r G S^, 
Jq 

/ CTh : e(v)rfx = F(v) for all v G Vh- 



For elements PS and ECQ4, the isoparametric bilinear interpolation is used for 
the displacement approximation, i.e. the displacement space Vh is chosen as 

Vh = {^r e V : ir = v\k o Fk e QiiKf for all KgTJ. (3.11) 

In other words, for v = (u, f )^ G Vh with nodal values v(Zj) = (wj, t>j)^ on K, 
4 



E 



1=1 



Ui 
Vi 



Vo + Vi^ + V2r] + Vu^f] 



(3.12) 



where 



/ Uo 




\ 




/ 1 


1 1 


1 \ 


/ Ui 


Vl \ 


Ui 


Vi 




1 


-1 


1 1 


-1 


U2 


V2 


U2 


V2 




~ 4 


-1 


-1 1 


1 


Us 


V3 


\ f/12 


Vl2 


I 




I 1 


-1 1 


-1 / 


\ U4 


f4 / 



We denote the symmetric stress tensor r 



For convenience 



Til ^12 
^■12 T22 

we abbreviate it to r = (th, T22, t^)'^- In [fT8l . the 5-parameters stress mode on 
K for the PS finite element takes the form 

/ 1 
1 




V 



4t \ 
73"^ \ 







r for/3^ := 



.AO 



VO 1 ^ry 



(3.13) 
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Then the corresponding stress space for the PS finite element is 

Sf^ := {r G S : f = t\k o Fk is of form flTT3|) for all K e Th} . 
In (301, the 5-parameters stress mode on K for element ECQ4 has the form 




62 
11 ' 



"1202 
"2 



^ ^ ^ \ 



t»12 
ai ' 



aibl2-ffli2fei 



6f £. 



202 C l_h2Ca^^ 6l„ 
62 ^ 02 ^ ai ' ai ' 62 



/3" for G 
(3.14) 



Then the corresponding stress space for the ECQ4 finite element is 



^EC 



{r G S : T = t\k o Fk is of form (llUl) for all K eTh] . 



Remark 3.4. T/ze stress mode ofECQ4 can be viewed as a modified version of PS 
mode with a perturbation term: 



1 "12 C 

b2^ 



012^2 



felfcl2 



ai2 , 



bi2 



01262-02612 t 
"2 

01612-01261 ^ 



ai2 
62 



2 

f'l , 



"2f \ 
62S \ 



01 

1 r/ 
1 4r/ e 



^ 62 01 ' 01 ' 62 ^ ' 

6?^ \ 



1 



hi 



+ 



/ — ,^ "12^12 01262-02612 ^ g g \ 

61612 012 

Oj ' Ol u,-^ 

612 ^ ai2 C 6i2f Oil 



'■-T] Slhz^mbi^ g g 

62^^ -t^'-^^ 0/ 



Remark 3.5. W/zen K & Th is a parallelogram, the stress mode of ECQ4 is re- 
duced to that of PS due to ai2 = &12 = 0. Thus, PS and ECQ4 are equivalent on 
parallelogram meshes. 



Define the bubble function space 

:= K e {L\n)f : v^(e,r/) = v^o^i^ e span{^'' -I , rf -if for all G Th]. 



Then for any \^ E B^, we have 



K 



(3.15) 
(3.16) 



with M^, ti^, f^, G R. 



Remark 3.6. It is easy to know (see l\26y} that for any K E Th, \u^ \ + \uri\ + \ v^\ + 



< 



Define the modified partial derivatives the modified divergence divv and 

the modified strain ^(v) respectively as follows [|34l : for K eTh, 

i -r dv , ^ X dy ,^ ^,dv dy ,^ ^,dv , dv , dv 



, ^ dv , „ , , ^ , dx , .dv dx , .dv dv dv 

~ I du dv 
divv A- = ^ + 

ox oy 

du if ^ _|_ 9v \ 

dx 2 9j/ 9x ' 



]_/&!£ _|_ 

2Vas/ aa;>' dy 



It is easy to verify that the PS stress mode satisfies the relation (see [|23l ) 

/ r : e(v'')dx = for all G S/^, (3.17) 
Jk 

or equivalently 

/ (r - To) : e(v^)rfx = for all v'' G S/^ 

for all T G S^"^, with tq the constant part of r, and that the ECQ4 stress mode 
satisfies the so-called energy-compatibility condition (see [[30l[35l ) 

[ T : £(v^)rfx = for all v'' G (3.18) 

JK 

for all T G Sf^. As a result, the stress spaces Sf^ can also be rewritten as 

Sf^ = {r G S : f^, = t,,|x o F,, G Pi{i,r]), [ r : £(v'')rfx = 0, 

Jk 

1 < « < J < 2 for all e B^, K e Tf,}, (3.19) 



Sf^ = {TGS:f,, = t,,|koF,, GPi(e,r/), / T:£(v^)dx = 0, 

1 < « < j < 2 for all E B^, K e Th}. (3.20) 

With the continuous isoparametric bilinear displacement approximation Vh given 
in (13.1 II) . the corresponding hybrid finite element schemes for PS and ECQ4 are ob- 
tained by respectively taking = and S/^ = Sf*^ in the discretized model 
(l3:9l ) (l3T0l) . 

Remark 3.7. Since the stress approximation of the hybrid elements is piecewise- 
independent, the stress parameters, f3'^ in A3.13\) or A3.14\) . can be eliminated at the 
element level. In this sense, the computational cost of the hybrid methods is almost 
the same as that of the isoparametric bilinear element. 
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3.3. Numerical performance of hybrid elements 

Three test problems are used to examine numerical performance of the hybrid el- 
ements PS/ECQ4. The former two are benchmark tests widely used in literature, 
e.g. |[l8l|T9l[30l[311[32l[35l, to test membrane elements while using coarse meshes, 
where no analytical forms of the exact solutions were given and numerical results 
were only computed at some special points. Here we give the explicit forms of the 
exact solutions and compute the stress error in L^-norm and the displacement error 
in iJ^-seminorm. For comparison, the standard 4-node displacement element, i.e. 
the isoparametric bilinear element (abbr. bilinear), is also computed with 5x5 
Gaussian quadrature. For elements PS and ECQ4, 2x2 Gaussian quadrature is 
exact in all the problems. 

E=1500 



2 



Figure 2: Cantilever beam 




regular meshes 





5 X 


1 











irregular meshes 



Figure 3: Finite element meshes 



Example 1. Beam bending test 

A plane stress beam modeled with different meshes is computed (Figure 2 and 
Figure 3), where the origin of the coordinates x, y is at the midpoint of the left end, 
the body force f = (0, 0)^, the surface traction g on r^r = {{x,y) G [0, 10] x 
[-1,1]: x = lOory = ±l}isgivenbyg|,,=io = (-2^2/, 0)^, g|j,=±i = (0, 0)^, 
and the exact solution is 



u 



—2xy 



-2Ey 




The displacement and stress results, ^"|^"''^^ and are listed respectively 

in Tables 1-2 with u = 0.25 and E = 1500. Though of the same first-order conver- 
gence rate in the displacement approximation, the hybrid elements results appear 
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much more accurate when compared with the biUnear element. Amazingly, the 
hybrid elements yield quite accurate stress results. 



Table 1 : The results of V'' in the plain stress beam test 

regular mesh irregular mesh 

method 5x1 10 x2 20 x4 40 x8 5x1 10 x2 20 x4 40 x8 

biUnear 0.3256 0.1106 0.03376 0.01165 0.5777 0.2668 0.09273 0.02881 

PS 0.07269 0.03635 0.01817 0.009087 0.1429 0.06303 0.03113 0.01552 

ECQ4 0.07269 0.03635 0.01817 0.009087 0.1313 0.06256 0.03107 0.01551 



Table 2: The results of ''"'n „° in the plain stress beam test 



method 




regular 


mesh 






irregular 


mesh 




5 x 1 


10 X 2 


20 X 4 


40 X 8 


5 X 1 


10 X 2 


20 X 4 


40 X 8 


biliear 


0.5062 


0.2951 


0.1545 


0.07826 


0.7242 


0.4854 


0.2809 


0.1481 


PS 














0.2663 


0.05559 


0.01134 


0.002551 


ECQ4 














0.1780 


0.03517 


0.007324 


0.001666 



Example 2. Poisson's ratio locking-free test 

A plane strain pure bending cantilever beam is used to test locking-free perfor- 
mance, with the same domain and meshes as in Figures 2 and 3. In this case, the 
body force f — (0, 0)^, the surface traction g on Fat = {{x, y) e [0, 10] x [—1, 1] : 
X = lOory = ±1} is given by g|x=io = {-'^Ey, 0)'^, g\y=±i = (0, 0)^, and the 
exact solution is 

_ f -2(1 - iy'^)xy \ _ f -2Ey \ 

'^~\{l-u^)x^ + u{l + u){y^-l) "^-y OJ- 

The numerical results with E = 1500 and different values of Poisson ratio u are 
listed in Tables 3-7. As we can see, the bilinear element deteriorates as z/ — )■ 0.5 or 
A — )■ oo, whereas the two hybrid elements give uniformly good results, with first 
order accuracy for the displacement approximation in if ^-seminorm and second 
order accuracy for the stress in L^-norm. 
Example 3. A new plane stress test 

In the latter two tests, the hybrid elements give quite accurate numerical results 
for the stress approximation. This is partially owing to the fact that the analytical 
stress solutions are linear polynomials in both cases. To verify this, we compute a 
new plane stress test with the same domain and meshes as in Figures 2 and 3. Here 
the body force has the form f = — (6|/^, 6x^)^, the surface traction g onF^r = 
{{x,y) : X ^ 10, -1 < y < 1} is given by g = (0, 2000 + 2y^Y , and the exact 
solution is 

i^+l. 4 4Nr f 2(a;M 

^-^iy^^) ' ^=[2{x' + y^) 
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Table 3: The results of for the bilinear element in the plane strain test 

Hi 1 



1/ 




regular 


mesh 






irregular 


mesh 




5 X 1 


10 X 2 


20 X 4 


40 X 8 


5 X 1 


10 X 2 


20 X 4 


40 X 8 


0.49 


0.9253 


0.7547 


0.4353 


0.1620 


0.8862 


0.7641 


0.5351 


0.2597 


0.499 


0.9921 


0.9690 


0.8866 


0.6619 


0.9515 


0.9241 


0.8530 


0.6978 


0.4999 


0.9992 


0.9968 


0.9874 


0.9514 


0.9615 


0.9567 


0.9446 


0.9067 


0.49999 


0.9999 


0.9997 


0.9987 


0.9949 


0.9626 


0.9606 


0.9591 


0.9540 



Table 4: The results of for PS in the plane strain test 

Hi 



V 




regular 


mesh 






irregular 


mesh 




5 X 1 


10 X 2 


20 X 4 


40 X 8 


5 X 1 


10 X 2 


20 X 4 


40 X 8 


0.49 


0.09759 


0.04879 


0.02440 


0.01220 


0.1557 


0.07342 


0.03649 


0.01822 


0.499 


0.09931 


0.04965 


0.02483 


0.01241 


0.1567 


0.07410 


0.03684 


0.01839 


0.4999 


0.09948 


0.04974 


0.02487 


0.01244 


0.1569 


0.07418 


0.03688 


0.01841 


0.49999 


0.09950 


0.04975 


0.02488 


0.01244 


0.1569 


0.07418 


0.03688 


0.01841 



Table 5: The results of ''"'h "^.^ for PS in the plane strain test 



V 




regular 


mesh 






irregular 


mesh 




5x1 


10 X 2 


20 X 4 


40 X 8 


5 X 1 


10 X 2 


20 X 4 


40 X 8 


0.49 














0.2286 


0.04566 


0.009326 


0.002094 


0.499 














0.2268 


0.0452 


0.009238 


0.002073 


0.4999 














0.2266 


0.04516 


0.009229 


0.002071 


0.49999 














0.2266 


0.04516 


0.009229 


0.002071 



Table 6: The results of '"iJ^^'''^ for ECQ4 in the plane strain test 



V 




regular 


mesh 






irregular 


mesh 




5 X 1 


10 X 2 


20 X 4 


40 X 8 


5 X 1 


10 X 2 


20 X 4 


40 X 8 


0.49 


0.09759 


0.04879 


0.02440 


0.01220 


0.1512 


0.07321 


0.03647 


0.01821 


0.499 


0.09931 


0.04965 


0.02483 


0.01241 


0.1526 


0.07392 


0.03682 


0.01839 


0.4999 


0.09948 


0.04974 


0.02487 


0.01244 


0.1527 


0.07399 


0.03686 


0.01841 


0.49999 


0.09950 


0.04975 


0.02488 


0.01244 


0.1569 


0.07418 


0.03688 


0.01841 



We only compute the the case of £■ = 1500, v = 0.25 for PS and ECQ4 and 
list the results in Tables 8-9. It is easy to see that the displacement accuracy in 
if^— seminorm, as well as the stress accuracy in L^-norm, is of order 1. 
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Table 7: 


The results of for ECQ4 in the plane strain test 








regular 


mesh 


irregular 


mesh 




V 


5x1 


10 X 2 


20 X 4 40 X 8 


5 x 1 10 X 2 20 X 4 


40 X 8 




A 
U 








0.1780 0.03456 0.007270 0.001661 


4QQ 


u 








0.1780 0.03455 0.007274 0.001662 


zLQQQ 


n 
u 








0.1780 0.03455 0.007275 0.001662 


n zLQQQQ 



u 








0.1780 0.03455 0.007275 0.001662 






Table 8: The error ^"i "'"^^ 

|u|i 


of Example 3 










regular 


mesh 


irregular 


mesh 




method 


10 X 2 


20 X 4 


40 X 8 80 X 16 


10 X 2 20 X 4 


40 X 8 


80 X 16 


PS 


0.1022 


0.05120 


0.02561 0.01281 


0.1815 0.08968 


0.04470 


0.02233 


ECQ4 


0.1022 


0.05120 


0.02561 0.01281 


0.1815 0.08968 


0.04470 


0.02233 


Table 9: The error ^^"'n of Example 3 






regular 


mesh 


irregular 


mesh 




method 


10 X 2 


20 X 4 


40 X 8 80 X 16 


10 X 2 20 X 4 


40 X 8 


80 X 16 


PS 


0.1022 


0.05120 


0.02561 0.01281 


0.1806 0.08590 


0.04239 


0.02113 


ECQ4 


0.1022 


0.05120 


0.02561 0.01281 


0.1850 0.09103 


0.04532 


0.02264 



4. Uniform a priori error estimates 
4.1. Error analysis for the PS finite element 

To derive uniform error estimates for the hybrid methods, according to the mixed 
method theory [|6ll3, we need the following two discrete versions of the stability 
conditions (Al) and (A2): 

(Al/t) Discrete Kemel-coercivity: For any r G := {r G S?, : /^^ r : e{w)dy^ = 
0, for all V G it holds 

Ikllo ^ a(T,T); 
(A2h) Discrete inf-sup condition: For any v G Vh, it holds 

|v|i < sup II II . 

O^reSh ||T||o 

Introduce the spaces 

Wh := {q G L^in) : q = q\K o Fk G Pi{k) for all K G T^, 
:= {qeWh-. q\K G Po{K) for all K G TJ. 
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To prove the stability condition (Alh) for the PS finite element, we need the fol- 
lowing lemma. 

Lemma 4.1. (134]) Let the partition satisfy the shape-regularity condition ( 13. il) . 
Assume that for any q G Wh, there exists some v G V/i with 

\\q\\l < / gdivvdx, < \\q\\l. (4.1) 



Then it holds 

g(divv + divv^)(ix 



g||o< sup -^""^ , ^ „ forallqeWh, (4.2) 



where the semi-norm | ■ on Vh + Bh is defined as \ ■ \ i,h '■= (Xlx 



I |2 



I 2 . 



Remark 4.1. Under the shape-regularity condition ([3.1\) , the following special 
property has been shown in ^34^: 

|v|i + \Ai,h < \^ + Ai,h < |v|i + Iv^li,;, for alive Vh,^r^ G B^. 

In view of this lemma, we have 

Theorem 4.1. Under the same conditions as in Lemma 4.1, the uniform discrete 
Kernel-coercivity condition (Al/jj holds for the PS finite element with (t^ = cr^^ . 

Proof Similar to the proof of Theorem 2.1, it suffices to show ||trT||o ^ ||'''^||o 
for any t E Zh. 

In fact, for t E Zh, for all v G V/i and for all E B^, it holds 



= I ''' ■ siv)dx. 
Jn 

= [ r: (£(v)+£(v^))rfx 
Jn 

= [ ( Jtrrl + T^) : (^(v) + e{A)d^ 
Jn ^ 

= / -trr (divv + divv^)rfx + / : {e{\) + e{A)dyi. 
Jn 2 Jn 

Thus, by Lemma 4. 1 , we get 

WtrrWo < sup ^ '—-^ '— 



-2Lx^:(e(v)+sV))dx 
^^P \ — : — M ^ If IIo- 



This completes the proof. □ 
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This theorem states that any quadrilateral mesh which is stable for the Stokes 
element Ql-PO is sufficient for (Alh). As we know, the only unstable case for Ql- 
PO is the checkerboard mode. Thereupon, any quadrilateral mesh which breaks the 
checkerboard mode is sufficient for the uniform stability (Alh). 

The latter part of this subsection is devoted to the proof of the discrete inf-sup 
condition (A2/j). It should be pointed out that in ||36l there has been a proof for this 
stability condition. However, we shall give a more simpler one here. 

From (13.121 ). for any v G V/i we have 



Jk 



dx 
ov 

dy 
du _|_ &u 
\ dy dx 



b2 + bl2C 





/ {Uib2 - U2h) + {Uih2 - Ui2bi)C + ([/1262 - U2bi2)r] \ 

{V2ai - Via2) + {Vi2ai - Viau)^ + (^2012 - Vi2a2)f? 

(C/2ai - Uia2) + (Uuai - Uiau)^ + {U2ai2 - C/i2a2)?? 
V +(Vi62 - V261) + (^1612 - Vi26i)e + (^1262 - 1^2^12)?? / 



-bi - bi2r] -61^ + ^2?? 

ai + aur] - a2'n 1 (4.3) 



-02 - ai2i ai + ai2r] ai^ - a2rj 



with/3- = {(3^ ■ ■ ■ , (3^,f := (f/i + ^Vi, U2 + !^V,, f/12 + '-^V,, V2 - ^V,, V,2 - 



Lemma 4.2. For any v G V/^ and K G Th, it holds 



\e{^)\\iK< 



mm 



(4.4) 



Proof. From (|4.3|) we have 

\H^)\\l,K = / ^(v) :e(v)dx 
JK 

[{{b2 + buOPl - (bi + bi2ri)Pl - {bii - b2^)(3lf + ((ai + ai2r,)Pl + (aiC 



+-(-(02 + ai2i)f3l + (ai + ai2ri)Pl + {aii - a2rj)^^ - (61 + bi2V)Pl " (^1^ 



- a2V)f3^f 
b2rjW,f 



< 



1 



l<j<5 



Lemma 4.3. For any r G and K G T^, it holds 



□ 



(4.5) 
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Proof. The form (13.131) indicates 



f T:Tdx= [ + vPl + )' + + + 

Jk J k L '^2 ^1 



•I? 



ai 02 



> - min JK(e,r/) (AO'- 



l<i<5 



□ 



Lemma 4.4. For anj v G Vh, there exists a E S^"^ such that for any K G Th, 



/ Tt, : £(v)(ix 
Jk 



\'^A\\k ^ lk(v)||o,if- 



(4.6) 



462 


-46i 








\ 













4ai 







-4a2 


4 02 J2 
'3-63- 


4ai 

3 ai 






4 Jo 
3 ai 
4 a2 Jo 

3~br 


-Ahi 

_AbiJi 
3 






4f)iJo 
3 

4 Jo 

3 62 / 





Proof. We follow the same line as in the proof of [Lemma 4.4, [fTOl l. 
For T G and v G Vh, from (I3TT31) and (|43]) it holds 



/ r : eiv)dx = (jS^f 
Jk 



By mean value theorem, there exists a point (^O) ^70) ^ [—1) 1]^ such that 

\\T\\lK = M^o,VoWf^(3^ 
with D = diag (4, 4, 8, |[1 + 2(^)2 + |[1 + 2(f )2 + (f^^ 

Denote D := diag (l, 1, 1, ^i^^^' 4+2^^+4 ^ 

( &2 









(4.7) 



-61 






















ai 
2 

_A 




Jo 


2 

feiJi 




bl Jo 


ai 




ai 
a2 Jo 


ai 




~«r 

Jo 
b2 



and take 



/ 1 



with 



1 Jry e 

Vo 1 f^y 



/3 



JK{^o,m) JK{^o,r]o) 
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•DA/3^ 



(4.8) 



we then obtain 

/ : e{v)dx = ||r„||^_^. (4.9) 
Jk 

On the other hand, (14^81) yields 
with 

/ 21 JL 2fei 

Jo 0.1 Jo Jfi 

a2 O1O2 262 



A 















ajb2 


-61 6i 





















Jo dlJo Jo 

aia2(fe2 Ji-fel J2) 2bife^ Ji _ a2bfb2Ji _ a2fef J2 2(ai&| Ji-a2bi J2) ajb2 
Jo^ J,f aiJ3 aija J» 

0^0 

aia2(-a2 Ji+ai J2) — 2a26ife2Ji 1 cq^^Ji , a2b\J2 2aia2(—b2J 

. P F o, 7^ W 

and D-i = diag (^1, 1, 1, "i+^y+^i ^ 4+2a^^y ^^^^ relation, together with 
Lemma[33 ^1} and dMl), imply 

l<i<5 1<«<5 

Combining this inequality with Lemmas [4 . 2H4 . 3 1 and (13.81 ). we arrive at 

ll_ Il2 > ll,/,,\l|2 



This inequality, together (14.91 ). shows the conclusion. □ 

Theorem 4.2. Le? partition Th satisfy the shape-regularity condition ( li.iD . T/zen 
uniform discrete inf-sup condition (A2h) holds with T,h = Sf'^. 

Proof. From Lemma l4!4l for any v G Vh, there exists r^, G Sf"^ such that (|4.6b 
holds. Then it holds 

||t^||o|v|i < j^T^ : T^dx^ j^ei^v) : e{^^)dx 

I '^v - Tydy. < I Ty -. £(v)(ix, 



< 



K 



where in the first inequality the equivalence of the seminorm |e(-)|o and the norm 
II ■ 111 on the space V is used. Then the uniform stability inequality {A2fi) follows 
from 

|v|i < J — < sup ■'^ ,, , ^ for all v G Vh. 

□ 
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Combining Theorem 4.1 and Theorem 4.2, we immediately have the following 
uniform error estimates. 

Theorem 4.3. Let (cr, u) G S x V" the solution of the variational problem 
d2.il)d2.2D . Under the same condition as in Lemma 4.1, the discretization problem 
d3.9Pd3.i0l) admits a unique solution [cTh, u^) G S^'^ x Vh such that 

||cr - CT/illo + |u - Uftli < inf ||cr - r ||o + inf |u - v|i. (4.10) 

In addition, let ph = ~\tr(Th be the approximation of the pressure p = — (/i + 
X)divu = —^trcr, then it holds 

\\p-Ph\\o^ inf Ik - '^llo + inf |u - v|i. (4.11) 

Remark 4.2. Here we recall that denotes "< C "with C a positive constant 
independent of A and h. 

Remark 4.3. From the standard interpolation theory, the right side terms of A4.10\) 
can be further bounded from above by Ch{\\(T\\i + ||u||2). 

4.2. Error analysis for ECQ4 

Since the stress mode of ECQ4 is actually a modified version of PS's with a pertur- 
bation term (see Remark 3.4), the stability analysis for ECQ4 can be carried out by 
following a similar routine. However, due to the coupling of the constant term with 
higher order terms, we need to introduce the mesh condition proposed by Shi [26J 
(Figure 4): 

Condition (A) The distance dxidK = "^^/du + ^12) between the midpoints of the 
diagonals of K E Th (Figure 2) is of order o(/ix) uniformly for all elements K as 




Figure 4: The distance dx 

For the uniform discrete kemel-coercivity (Alh) we need the following lemma. 

Lemma 4.5. (S34\l) Let the partition Th satisfy AS.U . Then for any q G Wh and 

V G Vh, there exists G Bh such that 

/ (g - Hog) (divv + dW)dx = ||g - Uoq\\l, (4.12) 
Jn 
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\Al,<\\q-Uoq\\l + 2\v\l (4.13) 
where Hq : L'^iVt) — )■ Wh is defined by Hoqlx '■= \ -JK^Qdx. 

We immediately have the following result. 
Lemma 4.6. Let the partition satisfy t\3.1\) and Condition (A). Then it holds 

, XX,, ,, ^ fo g (divv + div/iV^)c/x 

(l-o(l))||g||o < sup -^"^ I ,| ^ forallqeWh, (4.14) 

veVh.vftgSh |v + V \i^h 

where o(l) means o{l) ^ Q as h ^ and div/j denotes piecewise divergence with 
respect to T^. 

Proof. For any q G Wh, we can write 

q\K o Fk = q^ + gf e + q2V- 

Then it is easy to know that Hoglii' = q^ . 

By Lemma 4.1, there exists \ E Vh such that (14.11) hold with q = HqQ- On the 
other hand, from Lemma 4.3 there exists satisfying (I4.12I)(I4.13I) . 

Since it holds the relations 

/ Uoq divv''(ix = 0, 
Jn 

and 

{q — Hog) divv^c/x = / (g — Hog) divv^rfx, 
Jn 

it follows from (021) . (031) . (Ol) that 

|v + v''|i,?,||g||o + /f^nog divv^rfx 

< ||g-nog||g + ||nog||g + /f^nogdivv''c/x 

< J^iq — Hog) (divv + divv^)(ix + Hog diwdx + Uoq diw'^dx. 
= J^q (divv + divv^)(ix. 



(4.15) 



For the second term in the first line of (14.151) . from (13.161) . Remark 3.6, (13.11) . Con- 
dition (A), and Remark 4. 1, we have 

I [ Ilogdivv^fixl = ^|(6i2(m5 - m^) + ai2K - t^c))|nog|A-| 



^ ^ M.^±^Qy^| + i^^i + i^^i + |^;^|)||nog||o,i^ 

< ^(l)l^ + v'|i,/^||nog||o,K < o(l)|v + v^|i,;,||g||o, 

Ken 

which, together with (14.151) . yields the desired result. □ 
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From Lemmas 4.4 and 4.5 we know that, under the assumptions in the lemmas, 
the inf- sup condition 



||,||„< sup (divv + divv>.,) ^^^^„^^,y^ ,,,,, 

veVh,^rbeBh |V + V li,;, 

holds when the mesh size h is small enough. 

Therefore, following the same routine as in the proof of Theorem 4.2, we arrive 
at the following result. 

Theorem 4.4. Under Condition (A) and the same conditions as in Lemma 4.1, the 
uniform discrete kernel-coercivity condition (Alh) holds for ECQ4 with S/^ = Sf*^ 
and sufficiently small mesh size h. 

Next we show the discrete inf-sup condition (A2h) holds for the ECQ4 finite 
element. Notice that Condition (A) states 

max{|ai2|, I612I} = o(/iK), max{| Ji|, | J2I} = o(/i^). (4.17) 
Recall the element geometric properties (I3.7I) - (I3.8I) . namely 

ai ^ 62 ~ hx, max{a2,bi} < O^Iik), Jq ~ h]^. (4.18) 

This allows us to view all the terms involving one of the factors a^jbu, Ji, J2 
as higher-order terms. In this sense, the ECQ4 stress mode (13.141) is actually a 
higher-order oscillation of the PS stress mode (13.131) (cf. Remark [34l) . Thus, under 
Condition (A) Lemmas SSllMl also hold for ECQ4 stress space Sf^. 

As a result, we have the following stability result for the ECQ4 finite element. 

Theorem 4.5. Let the partition satisfy the shape-regularity condition di.iD and 
Condition (A). Then the uniform discrete inf-sup condition ( A2h) holds with 11^ = 

^EC 

Combining Theorem 4.4 and Theorem 4.5, we immediately have the following 
uniform error estimates for the ECQ4 finite element: 

Theorem 4.6. Let (cr, u) G ^ V be the solution of the variational problem 
^2.m2.2^ . Under the same conditions as in Theorem 4.4, the discretization problem 
d3.9Dd3.i0l) admits a unique solution [cTh, Uh) G Sf" x Vh such that 

||cr - cr/i||o + |u - Uftli < inf ||cr - t||o + inf |u - v|i. (4.19) 

In addition, let ph = —\tr(Th be the approximation of the pressure p = — (yU + 
X)divu = —^tra, then it holds 

\\p-Ph\\o< inf ||cr - r||o + inf |u - v|i. (4.20) 
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5. Equivalent EAS schemes 

By following the basic idea of [|2ni22ll23l . this part is devoted to the equivalence 
between the hybrid stress finite element method and some enhanced strains finite 
element scheme. 

The equivalent enhanced strains method is based on the following modified Hu- 
Washizu functional: 

n(T,v,7,7'') = + T : (7 - £(v) - y )cix 

- (b g-\ds- f ■ vdx}, 

Jfi^ndK Jk 



where 



b{a,/3) = a : C/3dx = / {2fia : f3 + Xtratrf3)d:s., 
Jn Jq 

V G Vft is the compatible displacements given in (13.111) . e{v) = (Vv + V^v)/2 
denotes the strain caused by the displacement vector v, r G S/t is the unconstraint 
stress tensor with 

:= {7 e L\n;Rlp : 7,^. = 7.,|xoF^ G span{l,^,ri} for z,j = 1,2, K e n}, 

7 G Eft and 7'' G are the independent strain and enhanced strain tensors respec- 
tively with 

Uj: = U'ps ■■= {s{^r') : G B,} 
for the PS finite element, and 

for the ECQ4 finite element. 

The variational formulations of the above enhanced strains method read as: Find 
{(Th, Uh, Sh, el) ethxVhxthxU^ such that 

V{/ T :(£;,- - s^)rfx} = for all tgS^, (5.1) 



6(7, Sh) - / 7 : (^hdyi = for all 7 G T.^, (5.2) 



n 



/ (Th : £(v)cix = V'{ / / ■ vdx + (h g • vds} for all v G \4, (5.3) 

Jn ^ JK Jr^ndK 



I ah : j'dx = 

n 
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for all 7^ G Ul 



(5.4) 



We claim that the hybrid stress finite element scheme (I3.9N3.10I ) for PS and 
ECQ4 is equivalent to the scheme (I5.1I) - (I5.4I) in the sense that the stress and dis- 
placement solution, (cT/j, Uh), of the latter enhanced strains scheme, also satisfy the 
equations (IX9l)(l3J0l) . 

In fact, we decompose T^h as T^h = ® 0-'h\^h), where = for the 
PS finite element and E,, = Ef*^ for ECQ4. It is easy to see that the relation (|54l ) 
indicates ah G E/i. Thus (15.41) is just the same as (13.101) . 

On the other hand, by using the decomposition of S/i, the equation (15.11) leads 

to: 

V{/ r : (e^ - eK) - s^)rfx} = for all r G E^\E/„ (5.5) 
K 

r-.ish- eiuh) - 4)rfx} = for all r G E^. (5.6) 

K 

Since 2fieh + XtrShl — cTh E S/j, from (15.21 ) we get dh = 2fieh + XtrShl or 
^ 271 ['''^ ~ 2(^+A) ^^^^-*-]- Substitute this into (|5.6I) . we then get an equation as 

same as (13.91 ). Hence, the equivalence follows. 

Notice that one can solve from the equation (|5.5I) . 

Remark 5.1. As shown in ^23\ |32]/ . we also have two higher-order hybrid stress 
finite element schemes equivalent to the schemes of PS and ECQ4, respectively. 
More precisely, the higher-order schemes are given as: Find {&h, i^h, ^1) ^ ^/j ^ 
Vh X Bh such that 

a{&h, r) - / T : {e{uh) + eM(u^)) cix = for all r G S^, 
Jn 

I &h : (£(v) + £m(v')) rfx = F(v) /or fl// v G \4, G B^, 
Jn 

where Em = ^ for the PS case and em = £ for the ECQ4 case. The equivalence is 
in the sense that the solutions of the scheme ( l3.9l) -( IX70l) for PS and ECQ4 and of 
the above higher-order scheme satisfy 

a-h = <^h and u,, = u/j. 

In fact, due to the constraints ^3.17^ - ^.18^ , we can view the higher-order scheme 
as an unconstrained one derived from the constrained scheme ( |3.9P -( |X7^ . with 
G Bh being a Lagrange multiplier 

Remark 5.2. Notice that in the hybrid stress finite element scheme f |3.9P -( |X70l ), a 
term like is involved. Thus for non-linear problems where C is not a constant 
modulus tensor, it is not convenient to implement the hybrid finite element method, 
while for the the enhanced strains method, this is not a difficulty, since one does 
not need to compute C^^. However, owing to the equivalence shown above, the 
hybrid finite element technology with PS and ECQ4 is easily extended to non-linear 
problems. 
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6. Uniform a posteriori error estimates for hybrid 
methods 



6.1. A posteriori error analysis 

By following the same routine as in Bl [8l HI, one derives the computable upper 
bound 

5Z ll^x(f + diva;,)||2_;^ + ||C-V/,-£K)||2_^+ Yl f^E\\[cThnE]\\lE 

(6.1) 

for the error ||cr — (Th\\1 + |u — VLh\\ of the hybrid finite element methods. Here £q 
denotes the set of all interior edges of Th, En the set of all edges on the boundary 
Fat, the length of an edge E E £ := SqIJ £n, ti^e ii^^ unit normal along E, and 
[cT/jii^;] the jump of CThn on E, especially for E G £n, [c/in^] := cr/jii^; — g. 
We first define an operator A : Tj x V — > (S x V)' by 

< A(cr, u), (r, v) >:= a(cr, r) — / cr : £(v)(ix — / r : e{u)dx. 

Jn Jn 

for all cr, T G S and u, v G V. Then, from (Al), (A2) and Theorem 2.2 we 
immediately get 

Lemma 6.1. The operator A defined as above is bounded and bijective, and the 
operator norms of A and A^^ are independent of X and h. 

We need the following weak interpolation operator 

Lemma 6.2. Let the partition Th satisfy d3.il) . Then there exists an operator J : 
V Vh such that, for all v E V, 

||/,^i(v - Jv)||o + \\h-'/\^r - Jv)||o,£ < |v|i. (6.2) 

In light of this lemma, we have the following a posteriori error estimate for the 
hybrid finite element scheme (I3.9I )- (I3.10I) . 

Theorem 6.1. Let the partition Th satisfy l lS.U . Then it holds 

||cr - cr/i||o + |u - u/j|i < ?7;,. (6.3) 

Remark 6.1. Here we recall that "<" denotes "< C "with C a positive constant 
which is bounded as X ^ oo and is independent of h. 

Remark 6.2. In fact, the reliable error estimate in Theorem 6. 1 is efficient as well 
in a sense that the estimate 



^ \\hK{i+d\vcTh)\\\j^+YhE\\[(ThnE]\\l^E ^ \\(T-(Th\\l+\u-nh\\+osc{i,Ti 

(6.4) 

holds, where osc(J,Th)'^ := Yl II ^/^(f ~ ^h)\\o k fa^ '^h piecewise constant 
integral means f/j. This can be obtained by following similar arguments in /EH/- 
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Proof of Theorem 6.1. The desired result can be obtained by following the same 
routine as in in [4]. Here for completeness we give a proof. 
In fact, the stability of A in Lemma 6.1 ensures that 

< A(cr - cr,„u- u/j), (r, v) > 
||cr - cr/t||o + |u - u/j|i < sup — — . 

res.vey Irllo + |v|i 

With the relation cr = C^^£:(u) and the Galerkin orthogonality j^{cr — cTh) '■ 
e{J'v)dx = 0, this equals 

/j^(C"^(cr - (Th) - e{u - Uh)) : T(ix - /^(cr - cTh) : £:(v)(ix 

sup j~j j~j ! ! 

tge.vgv IrWo + |v|i 

Ue{uh) - C^'^o-h) : T(ix - JJcr - (Th) : e{v - Jw)dyi 
= sup — j- — j- J— . 

xgs,vgv IrWo + |v|i 

With Cauchy's inequality and integration by parts, plus Lemma 6.2, this is bounded 
from above by 

sup - 5" / (f + divert) : (v - Jv)dx + V [cruris] • (v - Jv)(is /|v|i + 



6.2. Numerical verification 

We compute two examples. Examples 2 and 3 in Section 3.3, to verify the reliability 
and efficiency of the a posteriori estimator 7]h defined in (|6.1I) . We list the results of 
the relative error e^, the relative a posteriori error rjr, and the ratio r^r/cr in Tables 
10-12 and Figure 5 with 



„ 112 , |„ „ |2\l/2 
CT/illo + |U — Uh\i) 



iMl + \<r ' ''''' (ll-ll§ + |u|?)^/^- 

The numerical results show that the a posteriori estimator r]h is reliable and efficient 
with the ratio r]r/ being close to 1 in Example 2 and being around 4 in Example 3. 
It should be pointed out that in Figure 5 the mesh-axis coordinates 2, 4, 8, 16 denote 
the respective meshes 10 x 2, 20 x 4, 40 x 8, 80 x 16. 



Vh 
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Table 10: Numerical results of the a posteriori error estimator for PS in Example 2 







regular 


mesh 


of 


Figure 3 


irregular 


mesh 


of 


Figure 3 


V 




10 X 2 


20 X 4 


40 X 8 


80 X 16 


10 X 2 


20 X 4 


40 X 8 


80 X 16 




r?r(e-4) 


4.3306 


2.1653 


1.0826 


0.5413 


500.41 


99.415 


21.676 


4.6915 


0.49 


er(e-4) 


3.5126 


1.7563 


0.8781 


0.4391 


452.18 


93.579 


21.203 


5.1370 






1.23 


1.23 


1.23 


1.23 


1.11 


1.06 


1.02 


0.91 




rir(e-4) 


4.3300 


2.1650 


1.0825 


0.5413 


496.40 


98.740 


21.586 


4.6974 


0.499 


er(e-4) 


3.5331 


1.7665 


0.8833 


0.4416 


447.56 


92.648 


20.981 


5.0817 




Vr/ 


1.23 


1.23 


1.23 


1.23 


1.11 


1.07 


1.03 


0.92 




ry,(e-4) 


4.3300 


2.1650 


1.0825 


0.5413 


496.00 


98.677 


21.585 


4.7100 


0.4999 


er(e-4) 


3.5352 


1.7676 


0.8838 


0.4419 


447.10 


92.555 


20.959 


5.0764 




f]r 1 £-r 


1.22 


1.22 


1.22 


1.22 


1.11 


1.07 


1.03 


0.93 




??r(e-4) 


4.3300 


2.1650 


1.0825 


0.5413 


495.96 


98.671 


21.585 


4.7117 


0.49999 


er(e-4) 


3.5354 


1.7677 


0.8839 


0.4419 


447.05 


92.546 


20.957 


5.0759 






1.22 


1.22 


1.22 


1.22 


1.11 


1.07 


1.03 


0.93 


Table 11 


: Numerical results of the a posteriori error estimator for ECQ4 


in Example 2 








regular 


mesh 


of 


Figure 3 


irregular 


mesh 


of 


Figure 3 


V 




10 X 2 


20 X 4 


40 X 8 


80 X 16 


10 X 2 


20 X 4 


40 X 8 


80 X 16 




??r(e-4) 


4.3306 


2.1653 


1.0826 


0.5413 


480.69 


86.785 


18.365 


4.0010 


0.49 


er(e-4) 


3.5126 


1.7563 


0.8781 


0.4391 


359.44 


75.927 


17.426 


4.2483 






1.23 


1.23 


1.23 


1.23 


1.34 


1.14 


1.05 


0.94 




rir(e-4) 


4.3300 


2.1650 


1.0825 


0.5413 


480.66 


86.998 


18.514 


4.0744 


0.499 


er(e-4) 


3.5331 


1.7665 


0.8833 


0.4416 


359.37 


75.971 


17.436 


4.2495 






1.23 


1.23 


1.23 


1.23 


1.34 


1.15 


1.06 


0.96 




rir(e-4) 


4.3300 


2.1650 


1.0825 


0.5413 


480.66 


87.025 


18.538 


4.0941 


0.4999 


er(e-4) 


3.5352 


1.7676 


0.8838 


0.4419 


359.37 


75.977 


17.437 


4.2500 






1.22 


1.22 


1.22 


1.22 


1.34 


1.15 


1.06 


0.96 




?7r(e-4) 


4.3300 


2.1650 


1.0825 


0.5413 


480.66 


87.027 


18.540 


4.0965 


0.49999 


er(e-4) 


3.5354 


1.7677 


0.8839 


0.4419 


359.37 


75.977 


17.437 


4.2501 






1.22 


1.22 


1.22 


1.22 


1.34 


1.15 


1.06 


0.97 
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Table 12: Numerical results of the a posteriori error estimator in Example 3 







regular 


mesh 


of 


Figure 3 


irregular 


mesh 


of 


Figure 3 


method 




10 X 2 


20 X 4 


40 X 8 


80 X 16 


10 X 2 


20 X 4 


40 X 8 


80 X 16 




Vr 


0.4260 


0.2152 


0.1081 


0.05420 


0.6232 


0.3137 


0.1579 


0.0793 


PS 




0.1022 


0.0512 


0.0256 


0.0128 


0.1806 


0.0859 


0.0424 


0.0211 




T]r / Cr 


4.17 


4.20 


4.22 


4.23 


3.45 


3.65 


3.72 


3.75 




Vr 


0.4260 


0.2152 


0.1081 


0.0542 


0.5938 


0.3154 


0.1610 


0.0812 


ECQ4 




0.1022 


0.0512 


0.0256 


0.0128 


0.1850 


0.0910 


0.0453 


0.0226 




Vr/ ir- 


4.17 


4.20 


4.22 


4.23 


3.21 


3.47 


3.55 


3.59 



ratio ratio 




(a) Example 2: v = 0.49 (b) Examle 3 

Figure 5: The ratio rjr/ for PS and ECQ4 
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